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ABSTRACT 


The application of digital computers for the simulation of physical 
systems has become widespread. This paper describes a program designed 
to measure the frequency response of the simulation model of a linear 


system. 
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I. INTRODUCTION 


The application of digital computers for the simulation of physical 
systems has become widespread. Computer programs have been developed 
that simulate mechanical, electrical, and electro-mechanical systems, 
circuits, chemical processes, biomedical problems, traffic control, and 
SOron, 

One of the most important groups of programs for the simulation of 
continuous systems is that which employs digital-analog simulators; 
i1.e., programs which simulate the elements and organization of the 
analog computer. 

The eee. the investigation presented in this paper is to 
develop a technique -- a digital computer program -- for measuring the 
frequency response of the simulation model of a system. While such a 
program should be applicable to any digital simulation language, the 
program description presented in this paper utilizes International 
Business Machines Company's DSL/360 Digital Simulation Language. 
DSL/360, a System/360 FORTRAN IV program for the digital simulation 
of continuous system dynamics, employs the building-block approach of 
digital-analog simulators while providing the power of logical and 
algebraic equation notation. 

DSL/360 provides a basic set of function blocks from which a 
physical system may be modeled: integrators, limiters, pulse generators, 
function generators, and so on. In addition, FORTRAN library functions 
and functions from the Scientific Subroutine Package (SSP) may be 
utilized. In the event none of these satisfies the user's needs, the 


user may provide his own function blocks or subroutines. 


A digital computer program to determine the steady-state response of 
a system to a sinusoidal input is a useful adjunct to a digital simulation 
language. With it, non-linearities in an ostensibly linear system may 
be detected and evaluated. 


Reference 1 describes in detail the use and operation of DSL/360. 


IIT. DISCUSSION 


This chapter contains a discussion of the response of a linear 
system to a sinusoidal excitation and the basic concepts utilized in 


the development of a digital computer program to determine such response. 


A, SYSTEM RESPONSE TO A SINUSOIDAL INPUT 

When a system is linear, its response to a sinusoidal input is a 
sine wave of the same frequency (the higher harmonics are negligible). 
The magnitude ratio and phase of the response depend on the forcing 
frequency but not on the input magnitude [Ref. 2]. The condition that 
magnitude ratio and phase be independent of the te. amplitude is a 
condition for the linearity of the system. 

The steady-state frequency response of a stable, linear system to 
a sinusoidal input can be determined analytically from the system 
transfer function [Ref. 3]. The response to an input A sin wt is 
given by 

y = A| P(jw)|sin(wt + 6) 

where | P( jw) | is the magnitude of P(jw), $ is the argument of P(jw), 
and the complex number P(jw) is determined from the system transfer 
function P(s) by replacing the s with jw. The system output has the 
same frequency as the input and can be obtained Beal ecainene the 
input by P(jw) and shifting the phase angle of the input by the argu- 
ment of P( jw). | P( jw) | and ® for all w constitute the system frequency 
response, where P(jw) is the gain of the system for sinusoidal inputs 


with frequency Ww. 


B. RESONANCE IN A SYSTEM AND STABILITY 

The maximum value of the magnitude of the closed-loop frequency 
response of a system is a measure of the stability of the system 
[Ref. 3]. The frequency at which this maximum occurs is the resonant 
frequency of the system. 

It often occurs in reality that a linear system's response will 
contain non-linearities under certain conditions. As an example, a 
second-order system may experience system gain non-linearity as a 
result of resonance if the excitation frequency is at or near the 
system's natural frequency. Such non-linearity and any resultant 
effect on system stability are of interest to a design engineer, for 
instance, since the performance of the system with which he is con- 
cerned may be vitally affected. Being able to utilize a computer to 
determine system non-linearity during the design stage may result in 


savings of both time and money. 


C. DEVELOPMENT OF THE PROGRAM 

A program to determine the frequency response of a system to a 
sinusoidal input must first determine attainment of the steady-state 
condition upon application of the excitation and then calculate values 


of the steady-state response magnitude and phase. 


D. ASSUMPTIONS 
Severaliassumptions are made regarding the system, its excitation, 
and its response. 
Assumption 1. The system is linear and stable. 
Assumption 2. The system may be described in DSL/360 format. 
Assumption 3. There are no initial conditions; i.e., all initial 


Conditions are Zero at time zero. 
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Assumption 4. The excitation is described by 
A sin(wt + @) 
where A is the peak amplitude, w is the angular frequency in radians 
per second, and 9 is the phase (input phase is zero throughout this 
paper). 
Assumption 5. The steady-state output of the system is periodic 
with the same frequency as the input and may be described by 
B sin(wt + $) 
where B is the peak amplitude, w is the angular frequency in radians 


per second, and © is the phase of the output with respect to the input. 


E. DETERMINATION OF THE STEADY-STATE CONDITION 

In any digital simulation language, the independent variable is 
time. The computational process is an iterative one, performed at 
intervals of time specified by the user. Attainment of the steady- 
state condition is determined by a mechanical, value-comparing process 
rather than by the solution of an equation, such as that discussed 
earlier. 

The steady-state condition is determined in the program discussed 
in this paper by permitting the system to respond to the sinusoidal 
input until certain steady-state criteria are met. Ten cycles of input 
are permitted to pass before attempting any steady-state testing, by- 
passing computations involving the more wildly fluctuating transients 
and, thereby, decreasing computation time. The figure "ten" is an 
arbitrary choice, based on the assumption that the output will have 
settled to a condition approaching the steady state after ten cycles 


of input. 


ae 


~ 


The steady-state determination is accomplished by comparing the 
value of the output magnitude at each iteration point on one cycle 
with the value at the same point on a later cycle. When these values 
are reasonably close, steady state is considered to have been attained. 
To reduce the possibility of chance satisfaction of the steady-state 
criteria while the output is still in the transient or settling phase, 
comparison of output magnitude values is made in the program with not 
one but two later cycles, the third and the tenth. Again, these figures 


are arbitrary choices. 


F, DETERMINATION OF OUTPUT PEAK MAGNITUDE AND PHASE 

Once the steady-state testing commences, the program simultaneously 
computes the peak magnitude of the output and its phase with respect to 
thew npitee 

The peak magnitude of the output is nothing more than the highest 
steady-state value computed for the output magnitude. 

The phase is computed from the times at which an input cycle and 
its resultant output cycle cross their respective zero-value points. 


The equations used in the program are: 


—_ 2m(t.-t ) 
R T 

360(t.-t ) 
sees 1. 36 
D T 


where @. is the phase in radians, %. is the phase in degrees, t, is the 
R P P il 


D 
time at which the input cycle crosses its zero-value point, ys is the 


time at which the resultant output cycle crosses its zero-value point, 


and T is the period of both the input and the output. 


IZ 


III. THE PROGRAM 


This chapter presents a description of the program evolved to 
determine the response of a linear system to a sinusoidal excitation. 
A sample program of a single-run job is shown in Figure 3.1. Reference 
to this program is made in the explanations in the paragraphs that 
follow. The explanations presuppose some familiarity with DSL/360 on 


the part of the reader. 


A. THE SYSTEM 

The system used in the sample program of Figure 3.1 is second-order, 
with a natural frequency of ten radians per second, but any linear 
system would have served as well. The system is described, from its 
closed-loop form, as having a forward-loop transfer function of 
100/s(s + 1) and a negative unity feedback-path transfer function. 

The system might also have been modeled in its open-loop form. 


The system modeled in this manner is shown in Appendix A. 


B. CONSTANTS, PARAMETERS , AND PROGRAM CONTROLS 
Constants, parameters, and program execution controls are shown in 
the sample program following the title and general system description 
(in DSL/360, asterisks in column 1 indicate a comment card). 
The one constant, PI, is self-explanatory. 
The parameters must be determined by the user and are: 
OMEGA - the input frequency, in radians per second 
A------ the input peak magnitude 


FGAIN - the feedback-path gain 


1s: 


CALL ENDRUN 
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FIGURE 3.1 


Y-STATE RESPONSE TO A 


SAMPLE PROGRAM FOR 
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The program execution controls must also be determined by the user 
and are; 


FINTIM ------ the maximum simulation value for the independent 
variable, time 


DELT -------- the simulation interval; the unit of time for the 
integration routine to accomplish an integration step 


INTEG RKSFX - the fixed integration routine with an integration 
interval equal to one-half of the simulation interval, DELT 


FINTIM must be expressed numerically in DSL/360 and, in this program, 
is set arbitrarily at 10,000 seconds to ensure sufficient time for the 
response to attain the steady-state condition. Once steady-state is 
attained and the output peak magnitude and phase are calculated, the 
run is automatically terminated by a CALL ENDRUN statement. 

In this program, the execution control DELT is directly related to 
OMEGA, the input frequency. Since DSL/360 requires that DELT also be 
expressed as a numerical value, the DSL/360 user must calculate the 
value of DELT from the equation 
Ze | 


x PERIOD = N x OMEGA seconds. 


le 


DELT = 


The first iteration calculation is made at time zero, and an iteration 
calculation is made at discrete intervals thereafter, the length of the 
intervals depending on the jen eeacion scree being used [Ref. 1]. 

Since the output phase is not known in advance, assurance that an 
iteration calculation will be made at or er se near the peak- 

value point of each output cycle is made by setting a high value for 

N. Since the basis for steady-state testing is the comparison of values 
at each iteration point on one output cycle with the values at the 


corresponding point on two other output cycles, N should be an integer. 


It has been determined empirically that the value of N should be an 


i. 


integer no less than about twenty; this figure provides a trade-off 
between acceptable accuracy and reasonable computation time. A program 
for generating a table of DELT versus OMEGA, for N equal to twenty, is 


shown in Appendix B. 


C. THE INITIAL REGION 
The INITIAL REGION encompasses calculations, input and out put 
operations, and initializations that must be made once only at the 
beginning of a run. The values and equations shown in the sample 
program apply to all systems and need not be changed by the user. 
All variables shown in the INITIAL REGION in the sample program 


are either self-explanatory or explained more appropriately elsewhere 


in this chapter. 


D. THE DERIVATIVE REGION 

The DERIVATIVE REGION encompasses those calculations involving 
integration and derivatives of the state variables being integrated. 
The basic interval in the independent variable, time, for each pass 
through this region is the calculation interval. In this program, the 
interval is determined by INTEG RKSFX and is one-half the simulation 
interval, DELT. 

In the program, the DERIVATIVE REGION contains equations describing 
the input and system transfer functions and input-output relationships. 
While the DSL/360 function block TRNFR is used in the sample program, 
representing the system in either its closed-loop or open-loop form, 
the system could also have been represented by a set of ordinary dif- 
ferential equations [Ref. 1]. The use of ordinary differential 


equations does not alter the remainder of the program in any way. 
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STORAG, a DSL/360 translator command, and TABLE, a DSL/360 data 
input, are both associated with the use of the function block TRNFR, 


SINE() is the DSL/360 function block that models the input. 


E. THE DYNAMIC REGION 

The DYNAMIC REGION encompasses all time-dependent algebraic calcu- 
lations, other than derivative calculations for the integration routines 
in the DERIVATIVE REGION, plus all other operations which must be per- 
formed at each discrete value of time; i.e., all those operations which 
must be performed each iteration. It is in this region, in the sample 
program, that steady-state testing and response magnitude and phase 
calculations are accomplished. 

The variables in the DYNAMIC REGION of the sample program are: 


OUT ---- the value of the output magnitude at the point on the output 
cycle at which the current iteration calculations are being made 


DELOUT - the value of OUT, to be saved for comparison with the value 
of OUT at the same point on the third cycle following 


DILOUT - the value of OUT, to be saved for comparison with the value 
of OUT at the same point on the tenth cycle following 


EPS1 --- the allowable difference between DELOUT and OUT (the differ- 
ence between the current value of the output magnitude and the value at 
the same point on the third cycle preceding) within which DELOUT and OUT 
are considered to be equal 


MAXOUT - the maximum value of OUT; the peak value of the output 
magnitude 


DEROUT - the value of the slope of the output at the iteration point 

EPS2 --- the allowable difference between DELOUT and OUT and between 
DILOUT and OUT within which DELOUT, DILOUT, and OUT are considered to 
be equal 


EPS3 --- the maximum allowable value within which the slope of the 
output, DEROUT, may be considered to be at the peak of an output cycle 


AXIN --- the time at which the input cycle crosses its zero-value 
point 


ley 


XIN ---- the value of AXIN saved until the input cycle next crosses 
its zero-value point 


AXOUT -- the time at which the output cycle crosses it zero-value 
point 


XOUT --- the value of AXOUT saved until the output cycle next crosses 
its zero-value point 


XTDIFF - the time elapsed between XIN and XOUT 
RPHASE - the value of the output phase in radians 
DPHASE - the value of the output phase in degrees 


ARG ---- the value of the argument for the calculation of the output 
slope, DEROUT 


ERG ---- the value of the argument for the calculation of EPS3 
EPS1 is relatively coarse and is the criterion to be met by the 
output magnitude before the program will permit steady-state testing 
during any given iteration. EPS2 is a refined criterion for determining 
whether or not the train of cycles is sufficiently similar for the out- 
put to be considered in the steady-state condition. The run is continued, 
however, even after the EPS2 criterion is satisfied, until the peak value 
of the output magnitude has been attained, determined, in part, by HPS3. 
The numbers 605 and 2005 appearing in the arguments of the equations 
for DELOUT and DILOUT represent the maximum number of sampled values of 
OUT stored in meee intervals DELINT and DILINT, respectively. 
Reference 1 states that these ntmbers, which must be coded exp iter i 
as numerical, integer neat actin should be less than or equal to the 
delay interval divided by DELT. In practice, however, it was found 
that these numbers must be greater than the delay interval divided by 
DELT. For an N of twenty, then, in the case of DELOUT, this number is 


DELINT _ _3 x PERIOD 


DELT ~ 1/20 x PERTOD ~ °°” 
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For DILOUT, the number is 


DILINT 10 x PERIOD 
DELT 1/20 x PERIOD 


= 200. 

It can be seen, therefore, that the sample program will handle a pro- 
gram for any DELT where N is 204 or pean. If more than this number of 
iterations per cycle is desired, the integer numbers to replace the 
numbers 605 and 2005 may be determined from the above equations. 

ERG is the argument, in radians, for the slope of the output at 
the point on the cycle most remote from the maximum point which can 
possibly be obtained for the DELT specified. Ideally, an iteration 
calculation will occur at the maximum, and the slope will be zero. 
However, because the process is an iterative one and the phase is 
shifted by an indeterminable amount, an iteration calculation may occur 
as far away, in time, as one-half of DELT; i.e., the maximum point may 
be exactly half-way between two iteration calculation points. To 


ensure that the slope criterion, EPS3, will be no greater than the slope 


at this point, ERG is determined as follows: 


: : : 2 epi 
Ninn OF radians/iteration ace ee aS = OO 


Nr. of iterations/cycle 


7 2. x pa 
PERIOD/DELT 


a 2. Xo ple see 
PER LOD 


ERG 


pi _ 5 (2 x pi x a 
2 2 PERI OD 


DELT 
Oe | 0.5 ” PERIOD 


il 


Phase is calculated in the program so as to produce the result 
always as a lag indication by XTDIFF, else both the true phase angle 
and its supplement would be calculated. For a phase angle from zero 


degrees to 180 degrees, the program is satisfactory. However, because 


i 


all attempts to limit phase computations only to times at which the 
input and output cycles crossed their respective zero-value points in 
an upward (increasing value) direction failed, phase lag angles from 
180 degrees to 360 degrees (phase lead) are calculated as supplements 
of the true phase angle. To determine whether the program has calcu- 
lated true phase or its supplement, the user may either Bats: the 
asymptotes of the phase plot on a Bode diagram or, if the system 
description does not readily lend itself to this method, have the pro- 


gram provide a time graph of a few cycles each of the input and the out- 


put, from which phase lag or lead may be determined visually. 


F., THE TERMINAL REGION 

Entry into the TERMINAL REGION is made at the termination of a run 
for purposes of performing input and output operations, testing 
terminating conditions, processing results, changing parameter values 
and requesting a rerun, or terminating the job. In the sample program, 
the values of the variables and parameters desired in the printout are 
caused to be written in the TERMINAL REGION. 

If it is desired to make several runs in one job, each run at a 
different value of OMEGA, say, it is within this region that new 
values for the parameter OMEGA and the simulation interval DELT may 
be specified. Appendix A contains an example of such a multi-run job. 

A method for programming a multi-run job with different values for 
OMEGA for each run, calculated by the program, and for which DELT need 


not be changed for each run is discussed in Chapter IV. 
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IV. CONCLUSIONS 


This chapter presents a discussion of the results obtained, several 
limitations of the program as evolved, and recommendations for exten- 


sions of investigation. 


A. RESULTS 

Overall, the results of testing several linear systems were satis- 
factory. The response data provided by the program were consistently 
accurate, signifying that the steady-state testing scheme worked 
satisfactorily and consistently avoided chance satisfaction of the 
Steady-state criteria during the settling phase. 

Programs and results for several of the systems tested are shown 


in Appendix A. 


B. LIMITATIONS 

There are several limitations on the use of the program evolved, 
some inherent in DSL/360 and some a function of the value of parameters 
selected. 

Limitation 1. DELT must be calculated by the user and expressed in 
the program numerically. This could be eliminated if DELT could be 
expressed as a variable function of OMEGA, but this requires an altera- 
tion of DSL/360, beyond the scope of this paper. 

Basically, a new value for DELT must be calculated for each value 
of OMEGA used. However, the number of calculations for values of DELT 
may be reduced in a multi-run job by starting with a high value for 


OMEGA and a corresponding value for DELT (equal to, say, one-twentieth 


rae 


of the period) and including in the program a scheme to halve the pre- 
ceding value of OMEGA each run. DELT need not be changed, then, since 
each time OMEGA is halved, the number of iterations per cycle is 
doubled, and the criteria for DELT are satisfied (N remains an integer 
no smaller than twenty, doubling each successive run; the program as 
shown in the sample will accept any value of N to a maximum of about 
204). Appendix A contains an example of this method. 

Limitation 2, Since the independent variable is time, only time 
plots may be obtained without modification of DSL/360. Such plots may 
be desired if the user wants to determine whether phase lags or leads. 

Limitation 3. Use of too small a value for OMEGA, as compared to 
the value for the natural frequency of the system, may result in no 
program output. This is caused by the fixed integration called for by 
INTEG RKSFX. If this card is removed from the deck, the integration 
interval automatically defaults to that called for by the variable- 
step integration routine RKS. However, DELT then varies with the 
integration interval [Ref. 1], and the conditions required for steady- 
state testing may not be met because integration points on different 
cycles may not correspond. It was attempted to begin a run with the 
variable-step routine RKS and, after some relatively short interval, 
to impose the fixed-step routine RKSFX, These attempts were unsuccess- 
ful, however, because no matter where the INTEG RKSFX card was located 

-- even within a NOSORT region [Ref. 1] -- it was employed from the 


outset. 


C. RECOMMENDATIONS FOR EXTENSIONS OF INVESTIGATION 
There are several recommendations for logical extensions of this 


investigation. 


Zz 


Recommendation 1. Determine a method for calculating phase that 
will not result in computing the supplement of the phase angle, rather 
than the true phase angle, whenever the phase is leading (or, as 
written in the program, whenever the phase lags from 180 degrees to 
360 degrees). 

Recommendation 2. Alter DSL/360 to permit the designation of DELT 
as a variable function of OMEGA rather than as a numerical constant. 

Recommendation 3. Alter DSL/360 to permit graphing a variable 
versus some other parameter than time; specifically, it would be 
desirable to graph both MAXOUT and the phase versus OMEGA. 

Recommendation 4. Apply the basic concept of this paper to the 


measurement of frequency response for a non-linear system. 
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APPENDIX A 


This appendix presents six examples of systems used to test the 
program evolved for determining the frequency response. For each 
example, the system description, the program, and the frequency response 


is given. 
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EXAMPLE 1 


The system used in this example is identical to the system used in 
the sample program of Figure 3.1. This example shows one method of 
obtaining the response to more than one frequency in one job. DELT 
must be calculated and specified within the program for each OMEGA. 
The program is shown in Figure A.l, and the frequency response is 
shown in Figure A.2. 

For this example, the simulation time varied from 3.65 seconds to 
8.11 seconds for runs for nine different frequencies, or an average of 


6.05 seconds per run. The total execution time was 58.49 seconds. 


= 
CALL ENDRUN 


e 
(ea) e 
nae | e 
--— - 
more) = = 
a+ e ag nm 
zy o>) (ea) | 
” e =. =) r 
eR _ LE: cG oOo 
=O - s — oe 
Ioe ° <t <I ree 
zea | st Zs bi? ie Cc oo 
i! B ood Qa Qa — e 
‘tJ te ~ Re YW ke 
Pat AE ~ OO ke Cc Wu UL 40), Ss 
wet } 2g 1 = coat — In 
NCH (4 =~ a «f{ «I 8 < 
eal te Cn > t- - eS. <a 
Ih |] Or =F Cy m ¢f Ta) —Y gay er 
Jee Cc Oe Oe | b= t <I e 
Oe = us ea t-D = ~ 20 > a 
— Zu cr) z= * wroCaieas COad «FC i= $O 0 LJ 
ue N WJ eo eG F YG alae Ne” = Ce We Q V/) 
i‘ oo uw Oo OD { Ca tLLl bee =O) aS TlmO (ake e <I 
OW © — - © FFEZE I ann Vo een CC) ee -  o o  t  BO a 
Zul ee Q ~ FF HH Ne e OOO ~NeaweCag ocor 6ocG a. 
Oun AO ‘Ss © D HH tmeyCiu KO W em ONO WwiIOWa e Q 
AWN tl bot QO Zar KURMOY Meise GrROdKeEFEF AQNQNOYIWE o 
AZ Ze rad » eed CUO e Sat Hil 2S 2 NN eQonr LW 
wt eo <_{ ai <{eUwNwwy’ LG enmm a D *€4A 22000 uulirtrnwnw o —” 
Cae Ty (egy (alee (SDamZ2cez Wi etf\k-t- r= Ob Re ee OO KOOL OO EOI <{ 
- OQ wWwaDn~ wc ew FLNOLCD Se xX DDD ey mm OY 2 DMD aie 
= Fue = Ob =H eCZ NOOK OYVaQOO0O0Ou; Zeer CO14+HCae a 
lh Sie Omar OZN ee ON | OLKK Tam ef DeF II] ec eo aw 
r- b— <t oJ) “oa tig <I emt Com 8 LL we ee LUI UH <TD ie 0- ODO KK Ee OO o 
Ao ™m HO mLIA DYE ef eno DYYOD Ifte{( > FCrO\Y MMH KH HOODS aa 
KA INN eo OA t to zZOODR~— e <iq 0) = ot & ~ OUI f m= of eCiaf <{ ae 
WMI xA«lO %¥* elND “OwtiwS © 2S JH J <— MIRrtTZ”UM!] eZaOkU'Y © © Or 
CiQU\ ee e 0D eo | eee eOmty Pe aDO Cwiliwdie) OF ROGONZO hee OQ Z-RO 
(LOc<Im ot SUM Rie Le OVWCOMNMESVZEZU ZtHterOG OO BCC OUHO er ORKNINES *OMNNWO«ds 
IIcsrTost Oniuxt eee OU MN UU AN Lao eee e e ON OR RRo ZOD ol NW NUN ibWw QO DeK tao 
WEOns Ne RT re eR tOO POU WWI ROM NY OMNR TERE ENR clk COCO TOR INIOZLNUT TNO en 


ZOU 9k KOC ZZON WH UHOUNLMANMNDD *#FAEN | OZDIDHOTPOZDIDA UW WZ il KR KRU NNOM GS 0 emi), 
= WMI ZU ULI ee OY) Led KI OOK SG We tem COR AT ORF MOON Mm 240 Det Ke qd eke OD 
StU WENA OY VION POLI Kae Ware ON~ es ewe KM OF eM Sew Oe eS OTT Owen ZO 


“ZU W WIR COAN MOF AA Ke ZY ONO CWRU TOOL WoOna er KE Kr ao WU ow OUST KM 
QAO CIJAOOWW Ww KK YO OR RW OW OD mae WOU LULU €f ed eK OY 0 Qe em mI (OU) 
a <2 < DD mh on be Se 

Lu b- SO Oe > lus WW Litem TON on) tp} 
ad Yi <i Wik _ ro ee < @ ae) e OQ =Z 
a ZzuxwZzree rod Ouo2z Cy 
= Ada Ze "Wy — <I > Wi 
ee HOA Q NteO* * ot mee 


04487989 


Ael 


FIGURE 
FRIGRAM FOR SYSTEM OF EXAMPLE 1 


26 


~~ 

wt oO N aa) ~ Oo uN 
oO Ww) oO on oO > Oo 
oO . ®) uw oO oO 8] oO 
Ye) © = Va) ~~ Ye) aa) 
“N oO tf Va) ol = ~ s 
oO al = aD Oo Pg AJ 
cm oa) fa) N AJ NI N 
© © eO eO eO © «© 
ee eo Of etme Ne Me ve 
DO DO HO BO HO HO MO 


x zOw ZOe 20x zOer ZOrzOerzOS 
qFOVSANOZAOZAOSAOZ®AOSl AGS 
BOW OW OWO OO OW OW OWM 


(CONTINUED) 


FIGURE A.1l 


2/ 


PHASE 


Oe 
BOGe 


50 0001 





€ - ], saTdwexyg jo waqsks 10z3 asuodsay 


¢ V danola 


AQCNLINOVW 


Aouenbalg 


(O4S/davu) 
AONANOAYA 





danLINOVW 


2&8 


EXAMPLE 2 


The system used in this example is identical to the system used 
in the sample program of Figure 3.1. This example shows a second 
method, a variation of that shown in Example 1, whereby the response 
to more than one frequency may be obtained in one job. For this method, 
DELT need not be recalculated for each OMEGA, which is halved each run, 
Since halving OMEGA results in doubling the number of iteration inter- 
vals per cycle for a given DELT. To ensure that the criteria for DELT 
(and N) are always met, the highest value for OMEGA for which a response 
is desired in the fob 1S specified in the PARAM card. 

The program is shown in Figure A.3. The frequency response, 
identical to that obtained for Example 1, is shown in Figure A.2. 

Since the number of iterations per cycle doubles each run, simu- 
lation time for this method is greater than for the method of Example 
1, where the number of iterations per cycle remains constant for all 
frequencies (so long as the same value for N is used). For a maximum 
allowable simulation time of ten minutes (a local ruling specified for 
the computer on which these examples were run), the responses for a 
maximum of only six frequencies were obtainable on any single job, an 
average of 100 seconds per run. The method of Example 1 is, therefore, 
in general less costly. In general, simulation time may vary consider- 
ably from one system to another and, for any given system, from one 


frequency to another. 
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PROGRAM FOR SYSTEM OF EXAMPLE 2 
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EXAMPLE 3 


This example, too, uses the system shown in the sample program of 
Figure 3.1 and in Examples 1 and 2. In this case, however, the system 
is modeled in its open-loop form, rather than the closed-loop form 
employed in the sample and Examples 1 and 2. The program is shown in 
Figure A.4, and its frequency response, identical to those obtained in 


Examples 1 and 2, is shown in Figure A.2. 
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PROGRAM FOR SYSTEM OF EXAMPLE 3 
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EXAMPLE 4 


The system used in this example is open-loop, with a transfer 
function of 10/s(st+1)(s+5). The program is shown in Figure A.5, and 


the system response is shown in Figure A.6. 
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FIGURE A.5 
PROGRAM FOR SYSTEM OF EXAMPLE 4 
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EXAMPLE 5 


The system used in this example is open-loop, with a transfer 
function of (st+l)/(st10). To meet DSL/360 criteria for system 
modeling [Ref. 1], the division was carried out, yielding 1 - 9/(st10). 


The system was then modeled from the diagram shown in Figure A./7. 





FIGURE A./ 


Block Diagram of System for Example 5 
The program is shown in Figure A.8, and the frequency response is 
shown in Figure A.9. 

The phase of the output in this system leads the input. One 
limitation to the program, discussed in Chapter IV, is that the pro- 
gram computes the supplement of the phase angle rather than the true 
value of the phase whenever the output leads the input, as in this 


example. True phase is used in the plot of Figure A.9. 
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EXAMPLE 6 


The system modeled in this example is a magnetic tape transport 
mechanism which employs a vacuum column for tape loop control. The 
equations representing the system have been linearized for use with 
the program. 

The input is an AC signal to a DC motor, and the output is the 
velocity of the tape before it passes over the capstan. 

A diagram of the system is shown in Figure A.10, and the program 
containing the system modeling equations is shown in Figure A.1l. 

The frequency response of the system is shown in Figure A.12. 

It will be noted that a resonance peak occurs at a frequency of 
approximately 20,000 radians per second. All attempts to determine 
frequency response at or near this frequency were unsuccessful because 
the steady-state criteria were not met within ten minutes, the maximum 
simulation time available locally on the computer on which these 


examples were run. 
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APPENDIX B 


A program for generating a table of DELT versus OMEGA for an N of 
twenty is shown in Figure B.1l. To generate a table for any other value 
of N, the user need only specify the desired value of N in the program. 
The program may also be modified to obtain values for DELT for intervals 
of OMEGA and a maximum value of OMEGA other than the one radian per 
second and the 10,000 radians per second, respectively, shown in the 


program. 


C GENERATION OF TABLE OF DELT VERSUS OMEGA FOR N=20 


WRITE (6,1) 
1 FORMAT(1X,'TABLE OF DELT VERSUS OMEGA FOR N=20',//,1X, 
*'OMEGA' ,9X, 'DELT',/) 
PI=3.1415927 
OMEGA=1. 
AN=20. 
DO 3 I=1,10000 
DELT=2 .*PI/ (AN*OMEGA) 
WRITE(6,2) OMEGA, DELT 
2 FORMAT(F7.0,E18.7) 
OMEGA=OMEGA+1. 
3 CONTINUE 
STOP 
END 


FIGURE B.1 


Program for Generating Table of DELT versus OMEGA 
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